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Rapidly cooling a liquid may result in a glass transition, creating an amorphous solid whose shear 
and bulk moduli are finite. Even when done with constant density, these resulting moduli depend 
strongly on the rate of cooling. Understanding this phenomenon calls for analyzing separately the 
"Born term" that exists also in perfectly ordered materials and the contributions of the "excess 
modes" that result from glassy disorder. We show that the Born term is very insensitive to the 
cooling rate, and all the variation in the shear modulus is due to the excess modes. We argue that 
this approach provides a quantitative understanding of the cooling rate dependence of a basic linear 
response coefficient, i.e. the shear modulus. 



Introduction: The appearance of new amorphous 
solids like bulk metallic glasses in contemporary technol- 
ogy brings about a pressing need to develop a theoretical 
understanding of the effect of the protocols of prepara- 
tion on the resulting properties of the obtained materials 
P, It was excellently expressed in 0] that "Since 
this so-called glass transition is essentially the falling- 
out of equilibrium of the system because the typical time 
scale of the experiment is exceeded by the typical time 
scale of the relaxation times of the system, the resulting 
glass can be expected to depend on the way the glass 
was produced, e.g., on the cooling rate of the sample or 
the particulars of the cooling schedule" . Thus for exam- 
ple whether a bulk metallic glass will tend to fail via a 
shear-banding instability depends on how it was prepared 
■ But given the inter-particle potential and even the 
density of states, can we provide a theoretical framework 
to predict this dependence? 

In this Letter we focus on the linear elastic moduli of 
the produced amorphous solids, and provide a theoretical 
framework to understand their dependence on the rate 
of cooling. Since it is known that changing the material 
density has a well understood effect on the elastic mod- 
uli [3] , we concentrate here on cooling protocols that keep 
the density constant Q. In previous work the common 
approach to explain the protocol dependence stressed the 
local motifs, be them icosahedra, tetrahedra etc. 
but these did not provide a quantitative understanding of 
the issues at stake. Rather, we will argue in this Letter 
that one needs to distinguish between two mechanically 
significant features, one related to the volume averaged 
Born term (and see below for a precise definition) that 
is determined by things like density, average number of 
bonds per particle, strength of interactions etc., and an- 
other, non-affine term, which is all about the degree of 
heterogeneity in the material. This approach is also dif- 
ferent from the traditional view of structure and rigidity 
where one tries to explain the latter in terms of the long 
range correlations in the former Q- an impossible goal for 
most glassy systems. Instead we say that what matters 
are average properties related to density and compress- 
ibility and the degree of mechanical heterogeneities in the 



material. 

Numerical simulations: To prepare quality data for 
the present discussion we have performed 2-dimensional 
Molecular Dynamics simulations on a binary system 
which is an excellent glass former and is known to have 
a quasi-crystalline ground state [1, Q . Each atom in the 
system is labeled as either "small" (S) or "large" (L) and 
all the particles interact via Lennard Jones (LJ) poten- 
tial. All distances \ri—rj \ are normalized by Xsl, the dis- 
tance at which the LJ potential between the two species 
becomes zero and the energy is normalized by esL which 
is the interaction energy between two species. Tempera- 
ture was measured in units of e^^ / where ks is Boltz- 
mann's constant. For detailed information on the model 
potential and its properties, we refer the reader to Ref 
^ . The number of particles in our simulations is vary- 
ing between 400 to 10000 at a number density n = 0.985 
with a particle ratio Nl/N^ = (1 + V5)/4. The mode 
coupling temperature Tmct for this system is known to 
reside close to 0.325. All particles have identical mass 

mo ~ 1 and time is normalized to t^ = esL^SL^ /"nio- 
For the sake of computational efficiency, the interaction 
potential is smoothly truncated to zero along with its first 
two derivatives at a cut-off distance rc = 2.5. To prepare 
the glasses, we first start from a well equilibrated liquid 
at a high temperature of T = 1.2 which is supercooled 
to r = 0.35 at a quenching rate of 3.4 x 10"'^. Sec- 
ondly, we then equilibrate these supercooled liquids for 
times greater than 20ri.ci, where is the time taken for 
the self intermediate scattering function to become 1% 
of its initial value. Lastly, following this equilibration, 
we quench these supercooled liquids deep into the glassy 
regime at a temperature of T = 0.01 at various quench 
rates, from instantaneous to infinitely slow. The interme- 
diate quench rates were 3.2 x 10~^ • • • 3.2 x 10~^ in jumps 
of one order of magnitude. The infinitely fast quench was 
achieved by a conjugate gradient energy minimization of 
the equilibrated liquid at T = 0.35. The infinitely slow 
rate was replaced by taking the quasi-crystalline ground 
state as the reference state. One should appreciate that 
the slowest quench rate required 0.21 billion MD steps 
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FIG. 1. (Color Online). Typical stress vs. strain curve ob- 
tained for AQS straining of one realization of a system of 
A'^= 10000 particles but for different rates of quench, as pre- 
sented in the inset. Note the significant change in the shear 
modulus (the slope at 7 = 0) and of the yield peak stress 
where the system yields to plastic flow. 
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FIG. 2. (Color Online). Shear modulus at 7 = as a func- 
tion of quench rate. (N = 4900) In blue symbols x we show 
the shear modulus /i itself. The Born contribution is given 
by the red round dots. The nonafflne contribution is given 
in terms of the green squares. Of course the shear modulus 
itself is the difference between the two other terms. 



which translated to 7 days of CPU time for the largest 
system. 

Once we have the quenched sohds we can strain them 
using an athermal quasi-static (AQS) protocol to exam- 
ine their stress vs. strain curves. In each step of this 
procedure the particle positions in the system are first 
changed by the afhne transformation 
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This transformation results in the system not being in 
mechanical equilibrium, and we therefore allow the sec- 
ond step, a non afhne transformation — >■ -I- Ui which 
annuls the forces between the particles, returning the sys- 
tem to mechanical equilibrium. One should understand 
that the non-affine transformation is a direct result of the 
amorphous nature of the material: in a regular lattice 
without defects the afhne step would leave the particles 
in mechanical equilibrium. The resulting data for some 
representative quench rates are shown in Fig. [TJ We ob- 
serve that both the shear modulus and the yield peak 
stress (where the system yields to plastic flow) decrease 
significantly when the quench rate is increased. In Fig. 
[2] we present the shear modulus as a function of quench 
rate (see blue "x" symbols"). This is the phenomenon 
that we want to clarify in a quantitative way. 

Theory: By definition the shear modulus is the second 
derivative of the energy of the system with respect to the 
strain 7, i.e. 
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In our process the full derivative with respect to 7 
needs to be elaborated, Physically it is computed keep- 



ing the net forces zero on all the particles, since we move 
from one mechanical equilibrium state to another. Thus 
the derivative contains two contribution, one the partial 
derivative with respect to 7 and the other, via the chain 
rule, the contribution due to the non-affine part of the 
transformation idj- 151: 



d d ^ d dui _ d ^ d dui 
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where the second equality follows from the form of the 
non-affine transformation where dr; = dui. Applying 
this rule to the definition of /i wc end up with the exact 
expression [l3- 15| 
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where the first term is the well known Born contribution 
which we denote below as ^b- The second term exists 
only due to the non-affine displacement Ui and it includes 
the Hessian matrix H and the non affine "force" H llSl: 



dridrj 



(5) 

Needless to say, before we compute the non-affine contri- 
bution in Eq. (|4]) we need to remove the two Goldstone 
modes with A = which are the result of translation 
symmetry. 

It is very important to stress at this point that the sep- 
aration between the Born term and the non-affine term is 
not an arbitrary one. The Born term is very insensitive 
to the quench rate in our example, and this is usually the 
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case: it is only sensitive to average properties like den- 
sity, average number of neighbors and interactions [l6| . 
In Fig. [5] we show the result of calculating the Born term 
for all our samples as a function of the quench rate, (see 
red dots in Fig. [2]), and there is only minor dependence. 
This is not the case for the non-affine term, whose di- 
rect calculation is also shown in the same figure in green 
squares. We see that this term changes significantly, tak- 
ing upon itself the full blame of the change in the shear 
modulus as a function of quench rate. The sum of the 
two terms agrees to very high accuracy with the direct 
measurement of the shear modulus from the slopes of the 
curves in Fig. [T]at 7 = 0. 

The density of States: as said, the Born term is 
almost independent of the quench rate, and its value is 
very close to that of the reference state which is the quasi- 
crystalline ground state. To understand the non-affine 
term we need to focus now on the density of states ^'(A) 
where Ai are the eigenvalues of the Hessian matrix. For a 
purely elastic piece of matter lacking of any disorder we 
know that the density of states is determined by the De- 
bye theory, and in terms of the eigenvalues of the Hessian 
matrix we expect a constant density 
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for a purely elastic medium. (6) 



For all our finite quench rates we have disorder in the 
resulting solid, and accordin gly we expect to see excess 
modes at small values of A jlSl . [Toj . These modes are 
sometime referred to as the Boson peak [2^. Their den- 
sity of states is shown in Fig. [3] as a function of the 
quench rate. We see very clearly that the density of ex- 
cess modes increases near A = as the quench rate is 
increased. For comparison we also show in the upper 
panel of Fig. [3] the constant density of states of a refer- 
ence elastic medium. Since the non-afiine term in Eq. |4] 
has the inverse of the Hessian, any increase in the density 
of states near A — >■ should have a strong effect on the 
shear modulus as is shown by the direct calculation. It 
is important to realize that the density of excess modes 
does not depend on the system size, and see for example 
the lower panel of Fig. |3] in which the density of states 
for the same cooling rate but for two different system 
sizes are superimposed. Thus one expects the same ex- 
cess modes also in the thermodynamic limit, and below 
we will see why the shear modulus computed in our small 
systems remains the same when N ^ 00. 

A relevant characteristic of these modes is their partic- 
ipation ratio, which provides a feeling as to how extended 
or localized the modes arc. Denoting the eigenvector as- 
sociated with an eigenvalue A^ as 'S'i, we use the following 
definition of the participation ratio: 
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FIG. 3. (Color Online). Upper panel: Density of eigenvalues 
of the Hessian matrix as a function of quench rate, normalized 
to unity. Data was collected in bins of size SX — 2.1. The 
Debye cutoff frequency Xd is shown by the arrow. In the 
inset we show the density of states for A < 100, a region that 
the non-afBne term in the shear modulus is very sensitive to. 
Lower panel: Density of eigenvalues of the Hessian matrix 
for the instantaneous quench at two different system sizes 
N = 4900 and N = 10000. Note that this has the highest 
D{X) at A ^ and that the density of states remains stable 
when the system size is changed. 
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where is the iih eigenvector projected on the jth 
particle. For fully extended modes this number is of 
0{N) whereas for localized modes it can be much smaller. 
In Fig. |3]we show the participation ratio of the modes ob- 
tained at a four different cooling rates. We see that the 
very first modes, including the Goldstone modes, have a 
participation ratio of the order of N. There is a little 
dip before a quasi-plateau. The modes associated with 
this dip are sometime referred to as the " Quasi-Localized 
Modes" (QLM). Lastly there are high eigenvalue modes 
which are very localized - these are the modes associated 
with Anderson localization. 

Which modes contribute?: the major question of 
interest for us at this point is which modes should be 
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FIG. 4. (Color Online). The participation ratio of the mode 
^'a with eigenvalue A at four different cooling rates as shown 
in the inset. The system size is A'' = 4900, leading to 9800 
different modes. Note that contrary to Fig. |3]here the A-axis 
was not binned. 
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FIG. 5. (Color Online). The convergence of the non-affine 
term in the shear modulus as a function of the upper cutoff 
in the sum over the modes. The calculation here is for a single 
realization. 



considered to quantitatively account for the non-afRne 
part of the shear modulus. It had been conjectured that 
maybe the QLM might contribute a dominant contribu- 
tion [l^ [21I, H^l ■ Next we provide a quantitative discus- 
sion of this issue. While the region near A —J- is impor- 
tant, it is not sufficient to account for the full non-affinc 
contribution. In fact all the modes up to the Debye cut- 
off are necessary to saturate the value of the non-affine 
term in the shear modulus. 

In Fig. [S] we present the computed non-affine term 
in the shear modulus where we use all the modes up 



to Af, excluding the Goldstone modes. In other words, 
we compute V^^ J^l^s '^1/ where at = "^k ■ 3. The 
dependence on A^ is approximately exponential with a 
typical scale of 130. Thus summing up to this range of 
A yields about 63% of the wanted quantity. To achieve 
an accuracy of 99% one needs to sum up to the Debye 
cutoff. All the excess modes are necessary to get the 
right answer. On the other hand we see that since the 
density of states does not change with the system size, 
the calculation of the shear modulus via this method will 
provide the correct shear modulus that pertains to the 
thermodynamic limit even though our systems arc very 
small. 

We note that connections between the Boson peak and 
softenin g of the amorphous material were mentioned be- 
fore 17-3 21, 24 1 . These connections did not produce 



however a quantitative statement of the type presented 
here. The obvious advantage of the present approach is 
its generality. The separation between the Born term 
and the non-affinc term is natural and robust, applying 
equally well to any example of amorphous solid. The 
conclusion of this study is that when one sees a strong 
dependence of shear modulus one should seek the expla- 
nation in the non-affinc rather than in the Born term. To 
remove any doubt that this important conclusion is not 
model-dependent we repeated the present study for the 
very different Kob-Andersen model [23| and the purely 
repulsive model [2^ and found again that the Born term 
is highly insensitive to the cooling rate. The non-affine 
term is determined by the low lying eigenvalues of the 
Hessian, but "low-lying" does not necessarily means the 
A — >■ range or even the QLM's. An accurately con- 
verged calculation of the shear modulus require all the 
excess modes up to the Debye cutoffi It is possible how- 
ever that higher order elastic coefficients may require a 
smaller range of eigenfunctions since they are more sin- 
gular in terms of the inverse of the Hessian. 
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